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ABSTRACT 



Context. Supersonic turbulence in the interstellar medium plays an important role in the formation of stars. The origin 
of this observed turbulence and its impact on the stellar initial mass function (IMF) still remains an open question. 
Aims. We investigate the influence of the turbulence forcing on the mass distributions of gravitationally unstable cores 
in simulations of isothermal supersonic turbulence. 

Methods. Data from two sets of non-selfgravitating hydrodynamic FLASH3 simulations with external stochastic forcing 
are analysed, each with static grid resolutions of 256^, 512'^ and 1024^ grid points. The first set applies solenoidal 
(divergence-free) forcing, while the second set uses purely compressive (curl-free) forcing to excite turbulent motions. 
From the resulting density field, we compute the mass distribution of gravitationally unstable cores by means of a 
clump-finding algorithm. Using the time-averaged probability density functions of the mass density, semi-analytic mass 
distributions are calculated from analytical theories. We apply stability criteria that are based on the Bonnor-Ebert 
mass resulting from the thermal pressure and from the sum of thermal and turbulent pressure. 

Results. Although there are uncertainties in the application of the clump-finding algorithm, we find systematic dif- 
ferences in the mass distributions obtained from solenoidal and compressive forcing. Compressive forcing produces a 
shallower slope in the high-mass power-law regime compared to solenoidal forcing. The mass distributions also de- 
pend on the Jeans length resulting from the choice of the mass in the computational box, which is freely scalable for 
non-selfgravitating isothermal turbulence. If the Jeans length corresponding to the density peaks is less than the grid 
cell size, the distributions obtained by clump-finding show a strong resolution dependence. Provided that all cores are 
numerically resolved and most cores are small compared to the length scale of the forcing, the normalised core mass 
distributions are close to the semi-analytic models. 

Conclusions. The driving mechanism of turbulence has a potential impact on the shape of the core mass function. 
Especially for the high-mass tails, the Hennebelle-Chabrier theory implies that the additional support due to turbulent 
pressure is important. 

Key words. Hydrodynamics — ISM: clouds — ISM: kinematics and dynamics — ISM: structure — Methods: numerical 
— Methods: statistical — Stars: formation — Turbulence 
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1. Introduction 

The observed supersonic turbulence in the interstel 
lar medium (ISM) plays an imp ortant role in the 



process of star formation (e . g.. iMac Low fc K 



2004 lElmegreen fc Scalol 120041: IScalo fc Elmegreen 
Ballesteros-Paredes etllf 120071: iMcKee fc Ostrikeil 



essenl 
2001 



20071) . 



The origin of this turbulence and the characteristics of 
different turbulence driving mechanisms are still a matter 
of debate. Due to the ability to counterbalance gravity 
globally and to provoke gravitational collapse locally 
turbulence is expected to have a strong impact on the 

shape of the stellar Initial Mass Function (IMF). The Joh nstone fc Ballvl 
high-mass end (m > 0.5 Mq) is typically observed to 
exhibit a powe r law of the form Af(m)dm oc m~°'drn 
(e.g., ISalneted 119551: iKrounal 120011: IChabried 120031: 



iLada fc Ladal 120031) . N{m) is the number of stars in the 
linear mass interval dm with ISaloeteil (|1955[ ) power-law 
index a — 2.35. Observations of dense cores embedded 
in star forming regions of turbulent giant molecular 
clouds show a similar but slightly shallower power-law 
distribution o f mass with an exponent a i n the range 
of 1.5 - 2.5 (jElmegreen fc Falgaronj Il996t iMotte et all 
I1998D . The origin of the IMF and the Salpeter power law 
for the high-mass tail are still under debate but there 



is a possible connecti on between the co r e mass function 
(CMF) and the IMF (iMotte et al.lll998t iTesti fc Sargend 
1998; Wilha ms et a ll 120001: iJ ohnstone et al.1 l200 0l. '20011: 

Nutter fc Ward-Thompsonl 



2007t lAlves. Lombardi. fc Ladal 1200^" see however, th e 
critical discussion by Clark. Klessen. fc Bonnell |2007[ ). 
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iPadoan fc Nordlundl (|2002D proposed a theoretical ex- 
planation for the CMF/IMF based on scaling properties 
of turbulence and the Jeans criterio n for gravitational 
instab ility. Combining the formalism of iPress fc Schechted 
()1974|) for cosmological structure formation with the notion 
of a turbulent Jeans m ass, a new analytic theory of the 
IMF was formulated bv iHennebelle fc Chabrieil (|2008D . 

Numerical simulations of turbulent molecular clouds re- 
portcd in the literature (for instance, Kritsuk et ah 2007t 



Hcnnc belle et al.ll2008HBaneriee et al.i r2009: Schmi dt et all 
20091 : iFederrath et al.l [2009a ) inject the turbulent energy 



via differen t meth ods, but apart from an early study 
by iKlessenI (|200l[ ). there has been no systematic anal- 
ysis of how different driving schemes affect the CMF. 
Possible physical mechanisms for maintaining ISM turbu- 
lence range from stellar feedback (supernovae, outflows, 
ionizing rad iation) to large-scale dynamical instabilities in 
the Galaxy (lElmeereen fc Scaloll2004l: iMac Low fc KlessenI 
120041: In orman fc Ferr aralll996^ , In this work we investigate 
the influence of compressibility of the turbulence forcing on 
the mass spectra of gravitationally unstable cores. 

We apply a clump- finding algorithm to the density fields 
from recent simulations of supersonic isothermal gas, using 
the FLASH3 code to solve the equations of hydrodynamics 
on static grids with resolutions of 256^, 512 ^ and 1024'^ grid 
points and pe riodic boundary conditions (jFederrath et al.l 
l2008bL I2009allbh . The turbulent energy is continuously in- 
jected by a specific force f on length scales corresponding 
to wavenumbers 1 < fc < 3, where k is normalised by 2tt/X 
for a box of size X. Each component of this force is mod- 
elled by an Ornstein-Uhlenbeck process (lEswaran &: Popd 
119881: ISchmidt et all l200fil 120091: IFederrath et al.l l2009a^ . 
We define the integral scale L by the mean wavelength 
of the forcing, i. e., L = X/2. The relative importance of 
solenoidal and compressive modes of the stochastic forcing 
is set by the parameter ( £ [0,1]. The simulations repre- 
sent two extreme cases: Purely solenoidal forcing (V • f = 0, 
^ = 1) on the one hand and purely compressive forc- 
ing (V X f = 0, C = 0) on the other driven to a RMS 
Mach number of Airms ~ 5.5. Turbulence becomes statis- 
tically stationary for t > 2 T, where T is the dynamical 
timescale. In order to take the intermittent nature of the 
density field into account we analyse each of the 81 density 
snapshots available in the time interval 2T < t < lOT to 
compute time-averaged core statistics. For a detailed de- 
scri ption of the simulation s etup and numerical techniques 
see IFederrath et al.l (|2009al) . Recent results of isothermal 
turbulence simulations indicate a strong dependence of 
the density field s tatistics on the compressibility of the 
stochastic forcing (IFederrath et all l2008bt I Schmidt et al.l 
I2009t IFederrath et al.l l2009al lbir Due to this effect it is a 
reasonable assumption that the resulting core mass distri- 
butions should show different properties for different tur- 
bulence driving mechanisms. Furthermore, the datasets en- 
able us to study in which way the effective resolution of the 
simulation affects the resulting mass distributions. 

This article is structured as follows. In the next sec- 
tio n, we will review th e anal ytica l theories for the CMF 



bvlPadoan fc Nordlundl (|2002D and IHennebelle fc Chabrieil 
(j2008[) . For the sake of clarity, we express the main results of 
these theories in a consistent notation. In Sect. [21 we discuss 
the parameters determining the global mass scale, which 
plays a key role in computations of the core mass func- 
tions. The clump-finding algorithm is briefly described in 



Sect, m and the resulting mass distributions are presented. 
Particular emphasis is put on the question of numerical 
convergence. The distributions obtained by clump-finding 
are compared to semi-analytic computations of the CMF 
in Sect. [SJ and the results are summarized and discussed in 
Section [51 



2. Analytical theories for the core mass function 

In this work, we follow IHennebelle fc Chabrieil (|2008D and 
write Af{M) = dN/dM, where N{M) is the cumulative 
number of cores with masses below M. The function Af{M) 
is also called the mass spectrum. Because of the wide range 
of different masses, it is useful to define the CMF by the 
number of cores per logarithmic mass interval. To that end, 
we define two dimensionless mass variables: 

— rn :— M/Mq is the mass relative to the solar mass. 

— M := M/Mj is the mass relative to the thermal Jeans 
mass at the mean density p. 

While m is useful to discuss the CMF in an observational 
context, M is the most natural choice for theoretical con- 
siderations. The linear mass distribution is related to the 
distributions with respect to logm or logM as followsf^ 



AA(A/) = (M0m)-i-^ = (M°M)-i- 



dlogm ' " ' dlogM 
Integration of over M yields the total number of cores. 



(1) 



tot - / AA(Af) dM 





f 

J — t 



diV 
dlogm 



d log m ~ 



dN 



(2) 



d log M 



^ d log M. 



2.1. Padoan-Nordlund theory 



iPadoan fc Nordltindl (|2002l ) based their analytical approach 
on the following assumptions: The core size is comparable 
to the thickness of isothermal postshock gas, the turbulent 
velocity fluctuations follow a power law, and the minimal 
mass that is unstable against gravitational collapse is given 
by the thermal Jeans mass. They showed that the number 
of cores per mass interval is given by an expression of the 
form 

I'M 

A^(M) cx M-(^+iM pdfj(M')dM', (3) 



where the parameter x depends on the scaling properties 
of turbulence. Integration of the probability density func- 
tion pdfj(M') of the Jeans mass for a given distribution 
of the mass density yields the fraction of cores with masses 
M' < M that are unstable according to the Jeans criterion, 
M' (X p~^l'^ . For cores of sufficiently high mass, this integral 
asymptotically approaches a constant. For this reason, the 
high-mass tail of M{M) obeys a power law with exponent 

Equation ^ can be written in terms of the probability 
density function pdf((5) of the logarithmic density fluctua- 
tion 5 — \og{p/ p). Since the Jeans stability criterion implies 



In this article, log denotes the natural logarithm. 
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logM = -^S, we have pdfj(M') = -2pdf(5)/M', and it 
follows that 

dTV ~ f°° 
— -^(xM-M pdf(<5)d(5. (4) 

dlogM J -2 log M 

Note that the upper bound on the core mass corresponds 
to a lower bound on the density fluctuation. The distribu- 
tion Q is invariant under a change of the global scales. 
Obviously, this cannot be the physical CMF if > Mtot, 
where Mtot := {'^L)^p is the total mass in the computa- 
tional domain (also see the discussion of the mass scale in 
Sect. [3]). For this reason, the Padoan-Nordlund theory has 
to be considered as asymptotic description for Mj <^ Mtot- 
The analytic calculation of the mass spectrum is straight- 
forward if the probability density function of the mass- 
density is lognormal, i. e.. 



pdf(<5) = 



1 



a/2 



: exp 



2a2 



(5) 



where 6 — — (7g/2. 

Originally, the power-law exponent x was obtained from 
the jump conditions for shocks in magnetohydrodynamic 
turbulence, and it was shown that the resulting slope of the 
hig h- mass tail is close to the Sa l peter slope x = a—1 = 1.35 
(see lPadoan fc Nordlunj2002D . iPadoan et all (|2007t ). here- 
after referred to as PN07, determined x also for isothermal 
hydrodynamic (HD) turbulence: 



a^HD 



2/3' 



(6) 



where /3 is the slope of the turbulence energy spectrum. 
Using data from a driven isothermal HD turbulence sim- 
ulation, PN07 calculated numerical mass distributions by 
means of clump-findin g for purely therma l support. The 
spectral index P = 1.9 (jKritsuk et al. I [20071 ) implies xnu = 
2.5 for the slope of the high-mass tail. This value recovers 
well the general trend of their numerical mass distributions 
(see Fig. 2 and Fig. 3 in PN07). However, they did not per- 
form time-averaging for their mass distributions, hence no 
estimate of the statistical significance of their results was 
provided. It was concluded that the high-mass tail is gen- 
erally too stiff for HD turbulence to be consistent with the 
observed CMF. 



Applying an approximation that excludes the largest 
cores, it is possible to generalise their approach to arbi- 
trary PDFs of the mass-density. For cores on a length scale 
R comparable to the integral scale L, the mass distribu- 
tion depends on the derivative of pdf (S) with respect to R, 
which, in turn, depends on the slope of the spectrum of the 
density fluctuations. It would be worthwhile to investigate 
the influence of the spectral properties of the density field 
on the CMF. Because of the difficulties in the numerical 
determination of the derivative of the J-PDF, however, we 
restrict our investigation to length scales R L, for which 
the general form of the mas s spectrum (equation (33) in 
iHennebelle fc Chabrierll2008D simplifies to 



7V(M) ~^p(AI 



dR ) dR 



exp(J)pdf(<5). (7) 



For a comparison with the Padoan-Nordlund theory, the 
case of purely thermal support is treated as follows. The 
core mass M = MjM is related to the logarithmic density 
fiuctuation S by the Jeans criterion for stability against 
gravitational collapse, and R is given by the thermal Jeans 
length. Hence, 



M = ^ = exp(-V2). 



(8) 



The dimensionless scale parameter is defined by i? :— i?/Aj, 

where Aj := ay^Cs/{Gp)^^'^ is the thermal Jeans length 
for the mean density (cs is the isothermal speed of sound, 
G the gravitational constant, and aj a geometry factor). 
Substituting the inverse mass-density relation ([5]) , it follows 
from equation ([7]) that 



dN 



d log M 



-^A/-3pdf(-21ogM). 



(9) 



We note that this expression for the CMF is fully deter- 
mined by the lower bound of the integral in equation 
and, generally, the high-mass tail does not obey a power 
law. 

If turbulence pressure contributes to the support against 
gravitational collapse, an implicit relation between the 
gravitationally unstable mass and the density fluctuation 
results: 



2.2. Hennebelle-Chabrier theory 



The theory of IHennebelle fc Chabrierl (|2008[) provides a 
general framework for mass distributions of cores on the ba- 
sis of the Press-S chechter statistical formalism . Other than 
in the theory bv iPadoan fc Nordlundl (|2002[ ). there is an 
inherent notion of a core scale R that is connected to its 
mass. Assuming a log- normal distribution of the gas den- 
sity, analytic formulae for the mass distribution M{M) can 
be derived. In the simplest case, M is given by the thermal 
Jeans mass for a given density. While the Padoan-Nordlund 
theory only uses thermal support, the Hennebelle-Chabrier 
theory considers a combination of thermal and turbulent 
support. On the one hand, turbulence is promoting star for- 
mation through the broadening of the density PDF, while, 
on other hand, it is also quenching star formation through 
the scale-dependent support it exerts on cores. 



M = R[l + MlR^"^ 



log 



1 + TW^E^" 



(10) 



The intensity of turbulence is specified by the Mach number 
A^, of turbulent velocity fluctuations on the length scale 



:= 



V3 



(11) 



where rj — {(3 — l)/2, and /? is the slope of the turbulence 
energy spectrum function in the inertial subrange. Since the 
turbulent pressure on the length scale Aj equals p{M*Cs)^, 
the parameter measures the relative signiflcance of tur- 
bulent vs. thermal support for cores of size ~ Aj. The turbu- 
lent Mach number on the length scale R is given by A^^i?**. 
Note that A4^,R^ is about unity if R is close to the sonic 
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length scale (|Schinidt et al.ll2009t iFederrath etall l2009a^ . 
The CMF including turbulent support is given by 



dN 



2p 1 + (1 - ri)MlR''^ 



dlogM M? R3[l + {2ri + 1)MIR^'^] 
X pdf \\og{M/R^) \ , 



(12) 



where R is numerically determined by the inversion of equa- 
tion (ITU| for a given value of M. 



3. Clump-finding and scaling 

For the calculation of the mass distributions we applied the 
implementation of the clump-finding algorithm by PN07. 
As a first step the algorithm divides the three-dimensional 
density field into k discrete density levels pi based upon the 
settings of the input parameters pmin and / = pi+i/ pi with 
i = 1 . . .k. The former parameter defines the minimum 
density level in units of the average density p while the 
latter defines the spacing between two adjacent levels. In 
the next step each density level is scanned for regions of 
connected cells with density values higher than the current 
density level. A connected region is counted as an object 

1 a core) if its mass exceed s the local Bonnor-Ebert mass 
Bonnoil [19561: lEbertlll955f ). The ratio of the Bonnor-Ebert 
mass to the solar mass is given by 



1.18c3 



MqG^/^P^/^ 

-2 



1000 cm-3 



-1/2 



T 
10 K 



3/2 



(13) 



where the gravitational constant G = 6.67 x 
10~^ cm'^ g~^ s~^, Cs is the isothermal sound speed, T 
is the temperature, and p, is the mean molecular weight 
of the gas. The number density n = p/(/xmH), where mn 
is the mass of the hydrogen atom, is averaged over the 
connected region. In the last step, each object which can 
be split into two or more objects is rejected. For a more 
detailed description of the clump-finding method, see 
PN07. 

For the clump- finding analysis, it is crucial to consider 
the degrees of freedom in the scaling of isothermal, non- 
self-gravitating turbulence. For a given temperature T and 
box size L, the hydrodynamical scales of the system are 
fixed, and the forcing magnitude determines the RMS Mach 
number A^rms of turbulence in the statistically stationary 
state. If self-gravity is not included, the flow properties are 
invariant with respect to the choice of the mean mass den- 
sity p. This is palpable if the compressible Euler equations 
are written in ter ms of the logarithmic de nsity 5 — hi(p/ p) 
(see, for instance. IFederrath et al.l l2009al ). Therefore, the 
mass scale is arbitrary. 

On the other hand, if gravitationally unstable cores are 
identified via postprocessing, then the value of jtibe (|13p 
resulting from the chosen values of n and T fixes the mass 
scale. In this paper, we characterize the mass scale by the 
total number A'be of Bonnor-Ebert masses with respect to 



the mean density contained in the simulation box: 



n^BE(p) 



1000 cm-3 



3/2 



10 K 



1.782 pc 



(14) 



where mBE(p) is the Bonnor-Ebert mass for the mean den- 
sity p, and n = p/lpmn). The ratio of the thermal Jeans 
length with the geometry factor aj = 1.18 corresponding 
to the definition of the Bonnor-Ebert mass to the integral 
length is given by 

' (15) 



L 



1/3- 

BE 



The Hennebelle-Chabrier theory accounts for the depen- 
dence of the core statistics on the parameter A'be via the 
factor p/M^ ~ N^j^/{2Lf in the distributions ® and P^ . 
Apart from that, variations of the mass scale correspond to 
changes of the scale-dependent Mach number, 



V3 



BE 



8 



(16) 



where characterizes the ratio of turbulent pressure to 
thermal pressure on the length scale Aj (see egn. fTT]) . In or- 
der to include the influence of the turbulent pressure on the 
stability against gravitational collapse in the clump-finding 
tool, we modified the stability criterion by replacing the 
thermal speed of sound in the definition of the Bonnor- 
Ebert mass with an effective speed of sound (see, for in- 
stance, iBonazzola et al.lll987() : 



Coff 



1^2 

2 core 



{v). 



(17) 



The turbulent velocity dispersion acoveiv) is computed from 
the mass- weighted RMS velocity fluctuation with respect to 
the centre-of-mass velocity for a particular core. 

The physical scales that were chosen by PN07 are 
n = lO'' cm-^, X = 2L = 6 pc and T = 10 K. In this case, 
iVBE « 1-2 X 10^ and X°j/L « 0.04. For the simulations with 
the highest resolution, we have A5 « 20.8A, where A is the 
size of the grid cells. As we shall see in the next Section, this 
introduces a serious resolution issue, because many cores 
extend only over a few cells for this parameter set. A fur- 
ther difficulty is that the above parameters are atypical 
for observed molecular cloud properties. According to the 
Larso n (198 ll) relation, n 3000 cm~3 (L/l pc)~^ (see also 
Fakarone et al.lll992l:lHever et al.ll2009D . This suggests that 
the assumed forcing scale L = X/2 = 3 pc is too large for 
molecular clouds of mean number density n = 10"* cm^^. 
In order to analyse the impact of the physical scales, we 
re-computed the core mass distributions for L — 0.6 pc, 
without changing the temperature and the mean number 
density. For these parameters, the number density is consis- 
tent with the Larson relation within a factor of two, which is 
reasonable given the scattering of observed molecular cloud 
properties (Falgarone et al. 1992). In this case, A^be — 10'^, 
and Aj/L « 0.2. We also note, however, that the validity of 
La rson- type relations for the mass density was questioned 
bv lBallesteros-Paredes fc Mac Low! (|2002D . 
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Fig. 1. The core mass distributions for the 512'^ solenoidal (a) and compressive (b) simulations as a function of the clump- 
finding algorithm parameter / (increasing from top curve to bottom curve) which sets the relative spacing between two 
adjacent density levels. Error bars contain the la temporal fluctuations and are only indicated for / = 1.04 for the sake 
of clarity. 



4. Core statistics 

We begin with a detailed analysis of the mass distribu- 
tions obtained for purely thermal support and the PN07 
mass scale corresponding to N-qe ~ 1.2 x 10^. In partic- 
ular, the tuning of clump-finding algorithm and influence 
of numerical resolution are considered in Section 14.11 and 
14.21 Finally, the mass distributions for a smaller number 
of Bonner-Ebert masses in the box, A^be ~ 1-0 x 10'^, as 
well as the influence of turbulent pressure are discussed in 
Section I 



4.1. Tuning of the clump-finding algorithm 

As described in section [3] the clump-finding routine needs 
two technical input parameters. The relative spacing / be- 
tween two adjacent density levels and the minimum density 
level pmin above which the datacube is divided into discrete 
levels. The infiuence of these parameters was already de- 
scribed in PN07. For the minimum density level we choose 
the mean density po of the simulation box since we observe 
no significant differences in the resulting core mass distri- 
butions for values of /Omin < Po which was also found in 
PN07. The influence of the level spacing / on the result- 
ing mass distributions is shown in Fig. [T] The slope of the 
high-mass range and the total number of cores is very sen- 
sitive to the chosen spacing of the density levels. For larger 
/-values the algorithm detects less but more massive cores. 
For decreasing spacing between density levels and therefore 
finer "scanning" of the datacube, the massive cores are split 
into a higher number of lower mass cores which leads to a 
steepening of the high-mass tail of the core mass distribu- 
tions. For / = 1.04, 1.08 and 1.16, this effect lies within 
the temporal fiuctuations of the mass distributions as in- 
dicated by the error bars in Fig. [T] PN07 also concluded 
that the mass distributions are basically converged for 
/ < 1.16. For the calculations in the present study, we chose 
a value of / = 1.04. Nevertheless, comparing our numeri- 
cal mass distributions to the observed IMF/CMF should 



be done with caution, especially, in the high-mass range. 
Note that the parameter dependency of core mass distribu- 
tions is a com mon problem f or a?/ clump- finding methods in 
the literature. ISmith et al .' (2008*) applied different clump- 
finding methods to data of SPH simulations of molecular 
clouds with the result that the definition of an individ- 
ual core and its mass depends str ongly on the method and 
the parameter settings (see also iKlessen fc Burkertl[2000l: 
Klcsscn' l200lfl . In the observers co mmunity t he clump- 
finding routines C LUMPFIND (.Williams et all 1 1994) and 
GAUSSCLUMPS (iStutzki fc Guestenlll990D ^ commonly 
used to decompose position-position-velocity information 
of molecular clouds into individual cores. These tech- 
niques suffer from the same parameter dependency, as was 
shown bv Schneider fc Brooks, (,2004) and more recently by 
iPineda et al.l (|2009D . 



4.2. Numerical resolution study for purely thermal support 

Fig. [5] indicates that the mass distributions depend very 
sensitively on the grid resolution of the simulation. For 
A'be = 1.2 X 10^, the high-mass tails steepen systemat- 
ically with increasing grid resolution for both solenoidal 
and compressive forcing, showing no sign of convergence 
in the sense of an asymptotic limit, even for the highest 
resolution of 1024'^ grid points. The mass distributions for 
A'be = 10'^ appear to be better converged, at least in the 
case of solenoidal forcing. The statistical properties of the 
core mass distributions are summarised in Table [1] The 
mean mass m shifts to lower values with increasing reso- 
lution, approaching a time averaged value of (to) « 1 Mq 
for the 1024"^ runs, while the width am of the distributions 
decreases. This is caused by the fragmentation of the cores 
found in the low resolution simulations into smaller and 
denser cores in the high resolution simulations which also 
satisfy the condition for gravitational instability. The sum 
of all core masses, mcorcs, is a few percent of the total box 
mass TOtot- 
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solenoidal i^=1.0 

N^=256^ 

N- = 51/'' 

N^=1024^ 




0.10 r< 



0.01 



compressive ("=0.0 

N^=256^ 

N^=5'?^ 

N^=1024^ 




(a) solenoidal forcing, A'^be = 1-2 x 10''', thermal support 



(b) compressive forcing, A^be = 1-2 x IC, thermal support 



solenoidal ("=1.0 

N^=256^ 

N^=.'i1 
N^=1024^ 




compressive ("=0.0 

N^=256^ 

N^=512^ 

N^=1024^ 




^ 0.10- 



(c) solenoidal forcing, A^be = 10"^, thermal support 



(d) compressive forcing, A^be = 10"^, thermal support 



Fig. 2. Core mass distributions for solenoidal (left) and compressive forcing (right) at numerical resolutions of 256'^ 
(dashed), 512'^ (dot-dashed) and 1024"^ (solid), normalised to the total number of cores, A'tot- Error bars indicate la 
temporal fluctuations. 

Table 1. Time averaged properties of the mass distributions for solenoidal and compressive forcing calculated "with the 
physical scalings applied in PN07. 



solenoidal forcing, = 1.0 



grid points 


(m) [MoJ 




(ntot) 


{rricorcs/mtot) | 


%] XLS 


(2;mle) 






256^ 


4.04 ±0.18 


7.00 ± 0.69 


1540 ± 57 


5.1 ±0.2 


2.1 ±0.5 


1.8 ±0.5 


9.8 ±4.9 


0.03 ± 0.07 


512^ 


1.83 ±0.04 


2.26 ± 0.10 


3918 ± 103 


5.9 ±0.2 


3.1 ±0.3 


2.7 ±0.7 


6.1 ± 2.6 


0.06 ± 0.09 


1024=^ 


1.08 ±0.02 


1.05 ± 0.04 


7082 ± 219 


6.3 ±0.2 


3.1 ±0.2 


3.2 ±0.7 


3.5 ±0.9 


0.21 ±0.26 








compressive forcing, (" 


= 0.0 








256^ 


9.7 ±8.0 


34.36 ± 29.67 


590 ± 57 


4.0 ±0.7 


1.0 ±0.3 


1.1 ±0.3 


9.7 ±8.0 


0.05 ± 0.07 


512^ 


2.41 ±0.19 


5.54 ± 1.48 


1859 ± 138 


3.7 ±0.3 


1.6 ±0.2 


1.6 ±0.3 


5.2 ± 2.7 


0.18 ±0.22 


1024^ 


0.98 ± 0.04 


1.44 ±0.21 


4985 ± 429 


4.0 ±0.3 


2.1 ± 0.2 


2.2 ±0.3 


2.7 ±0.8 


0.23 ±0.21 



An important property of the mass distributions is the 
slope X of the high- mass "wings of dA^/dlogrn, "which is 
related to the exponent a of the linear distribution via 
X = a — \. Assuming that there is a po'wer la"w range, "we 
applied error-"weighted least squares fits to the high-mass 
po-wer-la-w regime (256^: m > ISMq, 512^: m > 8Mq, 
1024'^: TO > 3Mq) in Fig. [21 The po"wer-la"w exponents 
obtained by this method are referred to as xls- For the 



highest resolution we find Ils = 3.1 ± 0.2 for solenoidal 
forcing and a shallo"wer slope of a;LS = 2.1 ± 0.2 for com- 
pressive forcing. The mass distributions of the solenoidal 
runs tend to sho"w a higher degree of convergence. The dif- 
ferences in the slopes bet"ween the three solenoidal simu- 
lations decrease "with increasing resolution. The slope a;LS 
obtained from the 512"^ simulation is steeper by a value of 
1.0 than the 256"^ "while the values of xls for the resolutions 
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of 512'^ and 1024^ are nearly identical in the solenoidal case. 
In contrast, the slopes of the mass distributions obtained 
from the pure compressive forcing increase by a value of 
0.6 (256^ ^ 5123) ^ yg^i^g of q 5 (5^22 10243) 

indicating slower convergence than in the solenoidal forcing 
case. 

We also used the method of maximum likelihood esti- 
mation (MLE) (for details se e Appendix inspired by the 
work of IClauset et al.l (|2007f) to determine the slope xmle 
of the core mass distributions and to check if the distribu- 
tion really follows a power law. This met hod av oids biases 
that can result from least-square fits (e. g. lBauka.2007 ) and 
provides an estimate of the mass value mmin above which 
the power law assumption is made. Furthermore, the prob- 
ability p that the power-law assumptio n holds can be cal- 
culated. FoUowing lciauset et all (|2007t ). we rule out power 
laws if p < 0.1. The results of the MLE method applied 
to the data of our core masses are summarised in Table 
[TJ For the highest resolution we find a power law slope 
of XMLE = 2.2 ± 0.3 for m > mmin = 2.7 ± 0.8 Mq in 
the case of compressive forcing and xmle = 3.2 ± 0.7 for 
m > TOmin — 3.5 ± 0.9 Mq iu the solenoidal forcing case. 
The errors are the la temporal fluctuations of these values. 
The slopes are consistent with the values obtained via least 
squares fits. At numerical resolutions of 256^ and 512^, the 
time averaged p- values are less than 0.1 for both solenoidal 
and compressive forcing. Thus, they are not consistent with 
power laws. However, at a numerical resolution of 1024'^, 
we find p = 0.21 ± 0.26 for solenoidal and p = 0.23 ± 0.21 
for compressive forcing. Note that the standard deviations 
of these p-values are very large, such that instantaneous 
CMFs can either be consistent with power laws or do not 
exhibit power-law behaviour at the high-mass end of the 
distribution. We emphasise that it is absolutely necessary 
to compute time-averaged quantities in order to obtain sta- 
tistically meaningful results. Note that the p- value can only 
rule out the power law hypothesis while other distributions 
could be a better fit even for p > 0.1. The possible incon- 
sistency of the high- mass CMFs with a power law m odel 
was also discussed bv iBallesteros-Paredes et al] (|2006[ ). 

For N-BE — 10'^, the scatter around the tails of the mass 
distributions is too large for the MLE method to be applica- 
ble. Therefore, we only list the power-law exponents x^s of 
the high-mass tails of the CMFs in Tabled Although being 
quite inaccurate, especially due to the low number statis- 
tics, this method can be used to describe the general trend 
of the CMFs. For the solenoidal simulation with turbulent 
support and Nbe = 10"^ we do not give a value for the 
power-law exponent because the error bars are simply too 
large. Decreasing the total number of Bonner-Ebert masses 
in the box from Nbe from 1.2 x 10^ to 10"^ leaves the fitted 
value of XLS relatively unchanged (taking the large error 
bars into account). The effect of adding turbulent support 
on the slope of the CMF will be discussed in Section FOl 

In order to further investigate the effects of numerical 
resolution on the core properties, we look at the size r of 
the cores. We define the approximate size of a core as 

'^^'^cc^nsA, (18) 

where ricdis is the number of grid cells contained in a core 
and A = X/N with N = 256, 512, 1024 depending on the 
numerical grid resolution. We will first consider the mass 
scale chosen by PN07, i. e., A^be ~ 1.2 x 10^ Then we will 



Table 2. Least-square estimates of the power-law exponent 
xls obtained from the 1024^ simulations for different values 
of A^BE with/without turbulent pressure included in the 
core stability criterion. 



solenoidal forcing, ( — 1.0 





A^BE 


turbulent support 




0.04 


1.2 X 10" 


no 


3.1 ±0.2 


0.04 


1.2 X 10" 


yes 


2.8 ±0.2 


0.2 


1.0 X 10^ 


no 


2.5 ±0.7 


0.2 


1.0 X 10^ 


yes 




compressive forcing, = 0.0 


0.04 


1.2 X lO'' 


no 


2.1 ±0.2 


0.04 


1.2 X 10" 


yes 


2.1 ±0.2 


0.2 


1.0 X 10^ 


no 


2.0 ±0.4 


0.2 


1.0 X 10^ 


yes 


1.2 ±0.3 



comment on the changes arising from a lower number of 
Jeans masses in the box (corresponding to A^be ~ xlO'^). 
The typical size of large cores should be ^ Aj, which is 
about 5.2A for A^ = 256 and 20.8A for A = 1024. The frag- 
mentation properties of the gas on length scales less than 
about lOA, which encompasses most of the range of core 
sizes, are affected by numerical smoothing. Consequently, 
the probability density functions (PDFs) of r, which are 
plotted in Fig. |3l change substantially with resolution. In 
particular, one can also see that the fractions of cores 
with r > Aj are relatively small, especially for the high- 
resolution data. At lower resolution, however, these frac- 
tions increase, i. e., more big cores are found relative to 
the number of smaller cores. Thus, the maxima of the mass 
distributions shown in Figure [2] are shifted towards higher 
masses with decreasing resolution. In addition, the discrep- 
ancies in the high-mass wings might be caused by a bias 
of the clump-finding algorithm to select cores that contain 
several gravitationally unstable cores. 

For a certain variance, cr^, of the density fluctuations, 
we expect that cores of size ^ cr~^/^Aj are most frequent, 
because the Jeans length changes with the inverse square 
root of the mass density. Fro m the values of a for sole noidal 
and compressive forcing (see iFederrath et al.|[2009al Table 
1), it follows that (cTcomp/CTsoin)^/^ = \/5.9/1.9 « 1.76. This 
is in agreement with the relative peak positions of the dis- 
tributions for A^ = 1024 in Fig. [31 Moreover, the mean size 
f of the cores for solenoidal forcing is larger by about the 
same factor than r for compressive forcing, as one can see 
from the values listed in Table [3] 

The minimal core size is roughly given by the peak den- 
sities pniax in the turbulent gas. The definition of the Jeans 
length implies Tmin ~ (p/Pmax)^/^A j. Since p^ax is roughly 
500p for solenoidal forcing and lO^p for compressive forc- 
ing (jFederrath et al.ll2009bh . it follows that ri„i„ ~ 0.9A 
and Tmin ~ 0.2A, respectively, for the highest resolution. 
Consequently, the smallest cores are marginally resolved 
in the 1024'^ simulation with solenoidal forcing, while they 
are definitely below the resolution limit for compressively 
driven turbulence. As an indicator, we calculated the frac- 
tions fr<2 of cores with r < 2 (see Table [3]). Indeed, signif- 
icant fractions fr<2 are obtained for all resolutions in the 
case of compressive forcing. In the solenoidal simulations, 
on the other hand, fr<2 drops from 15.8% (256^) to 0.037% 
in the 1024'^ simulation. These trends are also visible in Fig. 

m 
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0.4 r 



0.3 r 



solenoidol i^=1.0 : 
n5=256- : 

N5=1024^ ; 




compressive ^=0.0 

N^=256J 

N-'=512-' 

N^=10243 




(a) solenoidal forcing, A'^be = 1-2 x 10^, thermal support 



(b) compressive forcing, A^be = 1-2 x IC, thermal support 



Fig. 3. PDF of core size r defined in eqn. ([T5|) for solenoidal (a) and compressive | (b) | forcing for a numerical resolution 
of 256"^ (dashed), 512'^ (dot-dashed) and 1024"^ (solid) and PMOV scaling. Error bars indicate la temporal fluctuations 
of the PDFs. 



Table 3. Time averaged properties of the core size PDFs for solenoidal and compressive forcing. 



solenoidal forcing, = 1.0 



grid points 


iVBE 


turbulent support 


(r) [cells] 


((Jr) [cellsj 


{.fr<2) m 


256^ 


1.2 


X 


lO'' 


no 


3.72 ± 0.08 


1.74 ±0.05 


15.8 ±1.6 


256^ 


1.2 


X 


10^ 


yes 


3.70 ± 0.08 


1.77 ±0.05 


16.7 ± 1.7 


256^ 


1.0 


X 


10^ 


no 


10.5 ±0.5 


3.4 ±0.3 





256^ 


1.0 


X 


10^ 


yes 


12.8 ± 1.0 


4.8 ±0.7 





512^ 


1.2 


X 


10^ 


no 


5.78 ±0.12 


2.51 ±0.05 


1.4 ±0.5 


512^ 


1.2 


X 


10^ 


yes 


5.82 ±0.13 


2.60 ± 0.06 


1.6 ±0.5 


512^ 


1.0 


X 


10^ 


no 


18.4 ±0.8 


6.1 ±0.5 





512^ 


1.0 


X 


10^ 


yes 


26.3 ± 3.0 


10.7 ± 1.9 





1024^ 


1.2 


X 


10^ 


no 


9.37 ±0.22 


4.17 ±0.10 


0.037 ± 0.044 


1024^ 


1.2 


X 


10^ 


yes 


9.82 ±0.24 


4.56 ±0.12 


0.041 ±0.050 


1024^ 


1.0 


X 


10^ 


no 


34.5 ± 1.8 


11.7± 1.1 





1024^ 


1.0 


X 


10=* 


yes 


60.8 ± 7.6 


23.1 ±4.3 





compressive forcing, = 0.0 


256^ 


1.2 


X 


10- 


no 


2.98 ±0.15 


1.67 ±0.11 


34.8 ±3.6 


256^ 


1.2 


X 


10^ 


yes 


2.94 ±0.15 


1.68 ±0.11 


35.9 ±3.6 


256^ 


1.0 


X 


10^ 


no 


6.0 ±0.6 


3.3 ±0.5 


2.9 ± 2.9 


256^ 


1.0 


X 


10^ 


yes 


5.8 ±0.6 


3.4 ±0.6 


3.5 ±3.6 


512^ 


1.2 


X 


10^ 


no 


4.18 ±0.19 


2.28 ±0.11 


14.7 ±3.8 


512^ 


1.2 


X 


lO'^ 


yes 


4.10 ±0.19 


2.29 ±0.12 


15.9 ±2.8 


512^ 


1.0 


X 


10^ 


no 


9.0 ±0.8 


4.7 ±0.9 


0.04 ±0.18 


512^ 


1.0 


X 


10^ 


yes 


8.9 ±0.9 


5.0 ± 1.2 


0.039 ± 0.2 


1024^ 


1.2 


X 


10^ 


no 


6.07 ±0.37 


3.42 ± 0.23 


3.7 ± 1.4 


1024^ 


1.2 


X 


10^ 


yes 


5.95 ±0.38 


3.45 ± 0.25 


4.2 ± 1.6 


1024^ 


1.0 


X 


10^ 


no 


15.1 ± 1.9 


9.0 ± 1.7 


0.008 ± 0.07 


1024^ 


1.0 


X 


10^ 


yes 


16.0 ± 2.4 


11.0 ±2.6 


0.012 ±0.10 



Setting A^BE — 10"^, the estimate of the minimal core size 
from the peak density yields rmin « 4.6A for solenoidal forc- 
ing and 1.0 A for compressive forcing. In the solenodial case, 
fr<2 drops to zero for the highest resolution. Compared to 
the PN07 setting, the cores are in general larger (see Table[3] 
and the top panels of Fig. S]) . Accordingly, the resolution 
dependence of the mass distributions is less pronounced, 
particularly, in the case of solenoidal forcing. However, com- 
paring t he ma ss distributions without normalization in Fig. 
5(a) and 5(c)[ one can see that the total number of cores de- 



creases by two orders of magnitude if the mass scale is cho- 
sen such that A^BE = 10'^. As in the case Nbe ~ 1-2 x 10^, 
we find that the mean core size for solenoidal forcing is 
roughly twice as large as for compressive forcing. 



4.3. Influence of the turbulent pressure 

As explained in Section [3l we used the effective speed of 
sound according to eqn. (jl7p in the definition of the Bonnor- 
Ebert mass, eqn. (fT5)) . to compute the mass distributions 
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0.15 - 



solenoidal ("=1.0 

N'=256^ 

N-=512-' 




0J5 - 



tk 0.10 - 



compressive ("=0.0 

N^=256' 

N^=5-2^ 
N^=1024^ 




(a) solenoidal forcing, A''be ~ 10^, thermal support 



(b) compressive forcing, A^be = 10'^, thermal support 



0.10 - 



0.05 - 



0.00 [Jl 



solenoidal ("=1.0 

N^=256^ 

N^=512^ 

N-'=1024^ 




(c) solenoidal forcing, A^be = 10'^, turbulent support 



0.10 



0.05 




compressive <f=0.0 

N^=256^ 

N^=512^ 

N^=1024^ 



(d) compressive forcing, A^be ~ W^, turbulent support 



Fig. 4. PDFs of core size distribution as in Fig. |3] but for Nbe = 10'^ with/without turbulent pressure support (see 
Sect. O and E]). 



of cores with thermal and turbulent support. For brevity, 
we will use the term turbulent support, for which it is un- 
derstood that the instability criterion is based on the sum 
of thermal and turbulent pressure. 



Fig. 5(b) and 5(d) compares the CMFs for TVbe = 1.2 x 
10^ and Nbe — 10"^ with turbulent support for the highest 
grid resolution of 1024'^. While turbulent support has no 
significant impact on the CMF for Nbe ~ 1-2 x lO'^, it 
clearly changes the CMF for Nbe — 10"^. For solenoidal 
forcing, turbulent support leads to a large shift of the peak 
in the CMF, corresponding to the pronounce d cha nge of the 
core size distribution. As one can see in Fig. 5(d)[ however, 
the temporal fluctuations are larger than the time average. 
This is a consequence of the small total number of cores. In 
the compressive case with turbulent support, the peak of 
the CMF is relatively robust, but we observe a flattening 
of the high-mass tail from xls ~ 2.0 for the CMF without 
turbulent support to xls ~ 1.2 (see Tabled, which is close 
to the Salpeter value of 1.35. 

Table [3] gives a summary of the basic properties of the 
core size distribution for all possible combinations of grid 
resolution, turbulent/thermal support and Nbe- If turbu- 



lent support is included, fr<2 becomes negligible both for 
solenoidal and for compressive forcing if = 1024. The 
core size PDFs are plotted in Fig. H) The large error bars 
are the result of the much smaller number of cores found 
in each snapshot compared to A'be — 1-2 x 10^. Including 
turbulent support, the probability to find large cores in- 
creases substantially for solenoidal forcing, while the PDF 
for compressive forcing does not change appreciably. In the 
case of solenoidal forcing, the largest effect on the core size 
is observed for the highest resolution of 1024^ grid cells. 
Thus, the support of cores against gravitational collpase is 
greatly enhanced by the turbulent pressure on length scales 
that are sufficiently large compared to the grid scale. In 
Sect. 15.21 we show that the turbulent pressure exceeds the 
thermal pressure on length scales above 30 A, i. e., outside 
the range of length scales that are subject to significant 
numerical dissipation. 



5. Comparison to theoretical models 

In this Section, the core mass distributions computed 
with the clump-finding algorithm are compared to theo- 
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Fig. 5. CMFs of compressive (solid) and solenoidal (dot-dashed) forcing for different values of A'^be with/without turbu- 
lent support (see Sect. 12.21 and 13)) and a grid resolution of 1024'^. The least-square fits to the high-mass tails are shown 
as dashed lines. 



ret ical predictions of Padoa,n fc Nordlundl (|2002D as well 
as iHennebelle fc Chabrieij ()2008[ ). The influence of the 
width of the densit y PDF resulting from d ifferent forcing 
was investigated bviHennebelle fc Chabrieil ( ^2009 ). Because 
iFederrath etall (l2008b|, l2009aD showed that the density 
PDFs are not exactly log-normal, we evaluated equa- 
tions (0]), dnj, and ([T2|) for the numerical PDFs calculated 
from the simulation data. T he PDFs for the three n umeri- 
cal resolutions are plotted in IFederrath et al.l (|2009aL Fig. 4 
and 6) . The resulting CMFs are thus semi-analytic. This is 
particularly important for the case of compressively driven 
turbulence, for which the core mass distributions calculated 
from log-normal fits to the PDFs of the mass density devi- 
ate substantially from the semi-a nalytic distributions b ased 
on the numerical PDFs (see also lSchmidt et aI1l2009l) . For 
the comparison with the data from our clump-finding anal- 
ysis, we normalise the distributions with the total number 
of cores iVtot, which is calculated from eqn. ©. Allowing 
for an arbitrary geometr y factor a,^ in the definition o f the 
thermal Jeans mass (see IHennebelle fc Chabri"eill2008l) . we 
set aj = f .18 and identify Mj* with the Bonnor-Ebert mass 



'7^BE(p)A^0 (see eon. I13|) . Thereby, we have a criterion for 
gravitational instability on the basis of the thermal pressure 
that is consistent with the clump-finding algorithm. 



5.1. Purely thermal support 

Substitution of the time-ave raged PDFs obtained from 
1024^ simulations fPedcrrat h'et al.l l2008bl l2009aD into 
equations ^ and ([5]), yields the semi-analytic CMFs plot- 
ted in the top panels (a,b) of Fig.|6l These distributions are 
independent of the parameter A^be because of the normal- 
ization by the total number of cores, iVtot- As one can see, 
both the Padoan-Nordlund and the Hennebelle-Chabrier 
theory imply that the peak of the CMF is shifted towards 
lower masses for compressively driven turbulence. This is 
a direct consequence of the shape of the density PDFs. 
Since compressive forcing produces larger density peaks, 
there is a higher fraction of small Jeans masses compared to 
solenoidal forcing. This results in the formation of smaller 
cores, as noted in Sect. 14. 2|) . As a consequence, the high- 
mass tail of the CMF (jH]) is significantly less stiff in the 
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MjMj M/Mj 
(c) turbulent support, X^/L = 0.04, iVBE = 2 x 10^ (d) turbulent support, Aj/L = 0.2, Nbe = 10^ 

Fig. 6. Comparison of the semi-analytic mass distribution following from the Hennebelle-Chabrier (HC) theory and 
the Padoan-Nordlund (PN) theory with the corresponding distributions obtained via clump-finding (large dots) for two 
different choices of global mass scale. The thin dotted lines are the tangents to the mass distributions with the Salpeter 
slope X = 1.35. For an explanation of turbulent support, see Sect. 12.2] and 151 



compressive case (the asymptote for M ^ 1 is mainly given 
by the factor pdf(-21ogM). The slope of the CMF (g]), on 
the other hand, is determined by the factor AI~^, where 
X is defined by equation ([6]), in the limit of high masses. 
The turbulence energy spectra computed from the simula- 
tion data yield the power-law exponents (3 = 1.86 ±0.05 for 
sole noidal forcing and = 1 .94 ±0.05 for compressive forc- 
ing (jFederrath et al.|[2009a|) . Thus, the slopes of the high- 
mass tails predicted by the Padoan-Nordlund theory are 
a; ~ 2.3 for solenoidal forcing and x ~ 2.7 for compressive 
forcing. 

The results from the clump-finding analysis are shown 
as dots in Fig. [HI (a,b). Let us first consider the case 
A^/L = 0.04 (Nbe — 2 x 10^). For compressively driven 
turbulence, the clump-finding data match the theoretical 
predictions very poorly. As discussed in Section r4.2[ the size 
of the smallest cores is less than the grid resolution in this 
case. Since the numerically unresolved cores are missing in 
the core statistics, the distribution is biased towards higher 
masses in comparison to the semi-analytic distributions. In 
the case of solenoidal forcing, on the other hand, the small- 
est cores are at least marginally resolved (see Table U]), and 
the low-mass wing as well as the peak position following 



from the clump-finding analysis are reasonably close to the 
theoretical predictions. While the distribution obtained by 
clump-finding roughly agrees with the CMF following from 
the Padoan-Nordlund theory also for M ^ 1, the slope 
of a; ss 3.1 (see Table [5]) is significantly steeper than the 
theoretical value 2.3. In Sect 15.21 it is shown that this 
discrepancy can be resolved by including turbulent pres- 
sure. Remarkably, the mass distributions from the clump- 
finding analysis are matched quite well by distributions of 
the form © for A°/L = 0.2 {Nbe = 10^). There are small 
deviations of the high-mass tails, which are, however, well 
within the error bars (see Fig [5]). In comparison to the 
distribution predicted by the Padoan-Nordlund theory, 
large discrepancies become apparent. Since the cores are 
sufficiently resolved both for solenoidal and for compres- 
sive forcing for Aj/L = 0.2 (see Tabled]), the most likely 
explanation is that A^be is too small and, consequently, the 
high-mass cores are not within the asymptotic regime, for 
which the Padoan-Nordlund theory applies. 
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M/M°j 




-15 -10 -5 5 

(a) AS/L = 0.04, TVbe = 2 x lO'' 



10 -5 5 

(b) Xyi = 0.2, iVBE = 10=^ 



Fig. 7. Normalised gravitationally unstable mass as a function of the logarithmic density fluctuation for purely thermal 
support (thick dashed lines) and with additional turbulent support (solid lines), as defined by equations ([5]) and (fTU|) . 
For solenoidal forcing (see Table |4|), the mass corresponding to the scale for which thermal pressure equals turbulent 
pressure (A^,7?'' = 1) is indicated by the dot-dashed horizontal lines, and the dashed horizontal lines specify the mass 

~ 1/3 

for which the core scale equals the box size {R — N^'^ ). 



Table 4. Dependence of various parameters on the mass 
scale and the forcing. 



solenoidal forcing, A^rms = 5.3, 


V = 


0.43 




Ms 


Mr=2L 


0.2 1.0 X 10^ 4.6 1.53 


0.74 


52 


0.04 1.2 X 10^ 0.9 0.77 


3.7 


256 


compressive forcing, Mnns = 5.6 


, V = 


-- 0.47 


0.2 1.0 X 10^ 1.0 1.52 


0.82 


57 


0.04 1.2 X 10'' 0.2 0.72 


4.1 


283 



5.2. Turbulent support 



Since iPadoan fc Nordlundl ()2002|) considered only the ther- 
mal Jeans mass, we will concentrate on the Hennebelle- 
Chabrier theory for the case including both thermal and 
turbulent support. Eliminating R from equations ([TU)) by 
means of numerical root finding yields the mass-density 
relations plotted in Fig. [T] For comparison, also the re- 
lation ([S]) for purely thermal support is shown. One can 
see that low-density cores of size R greater than Aj can 
maintain a much higher mass against gravitational col- 
lapse if turbulent pressure is included in the Jeans crite- 
rion. The high-density asymptote, on the other hand, coin- 
cides with the thermally supported branch, because cores 
of sufficiently high density are associated with small length 
scales R, for which the turbulent velocity dispersion be- 
comes negligible compared to the speed of sound, i. e., 
M^R'' < 1. This implies M ^ R exp(-(5/2). For the 

length scale R^ ~ L{Mr,r,s/V?>y^''\ we have M^R^ = 1, 
i. e., the turbulent pressure of the gas equals its ther- 
mal pressure. Substituting the values of the RMS Mach 
number s and the scaling exponents from iFederrath et al.l 
(|2009a|) , the values of listed in Tableware obtained. 
It follows that i?s « 0.074L and 0.082L for solenoidal and 
compressive forcing, respectively. These scales are close to 
the sonic length scales As ~ 0.077L and 0.074L, which 
IFederrath et al.l (|2009al ) determined from the turbulence 



energy spectra. The corresponding dimensionless masses, 
Ms, are also listed in Table SQ 

An important limitation of calculations based on the 
approximation ([7]) is that R has to be small compared to 

~ 1/3 

the box size 2L, i. e., i? ^ -^be • "^^^ mass parameters 

~ 1/3 

Mii=2L corresponding to i? = iVg^ are indicated by the 
dashed horizontal lines in Fig. [T] The values of Mr=2L are 
listed in Table Sj for all parameter sets. As one can see, we 
have the constraints M < 50 and M < 250 for Aj/L = 0.2 
and 0.04, respectively. Since these limits are much higher 
than the peak positions of the mass distributions, there is 
a sufficient margin to study the high-mass wings. 

The PDF data for the density fluctuation 5 yield the 
mass distributions plotted in the bottom panels (c,d) of 
Fig. |6l Comparing the mass distributions with and without 
turbulent support, the high-mass tails are flatter in the for- 
mer case. As expected, the difference is more pronounced 
for Aj /L = 0.2, because of the higher contribution of turbu- 
lent pressure for large core masses (see Fig. [7]) . Remarkably, 
the tail of the mass distribution for compressively driven 
turbulence that is plotted in Fig.[6](d) is much less stiff than 
what PN07 reported for hydrodynamic turbulence and is in 
very good agreement with the Salpeter power law. The peak 
positions, on the other hand, are nearly unaffected, because 
the turbulent pressure on the length scales corresponding 
to the peaks of the mass distributions is small compared to 
the thermal pressure. 

Regarding the mass distributions obtained by clump- 
finding with turbulent support (large dots in the bottom 
panels of Fig jH]), similar trends as for the distributions 



^ These values implicitly depend on the integral scale L, which 
is set to half the box size (see Sect.[T|). More accurately, L can 
be determined from the velocity structure functions (A. Kritsuk 
2010, private communication). Indeed, this would entail a some- 
what closer match between the semi-analytic core mass distri- 
butions and the clump-finding data. Since the trends remain 
unaltered, however, we keep the simple definition of L in this 
article. 
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with purely thermal support can be seen. For Aj/L — 

0. 04 (A^BE ~ 1-2 X 10^), the clump-finding distribution is 
shifted towards higher masses for compressively driven tur- 
bulence because the smallest cores cannot be resolved (see 
Sect. 15. ip . The discrepancy is much smaller in the case of 
solenoidal forcing, for which the overall shape of the semi- 
analytic and clump-finding distributions agree quite well. 
However, the slopes of the tails following from the clump- 
finding data are steeper in both cases. Consequently, it ap- 
pears that the clump-finding algorithm underestimates the 
turbulent support of high- mass cores, provided that the the- 
oretical description of the CMF is correct. 

In the case X^/L — 0.2 (A^be — 10'^), the mass distribu- 
tions agree well for compressive forcing (except for a small 
shift), whereas the mass distribution obtained by clump- 
finding is markedly different from the semi-analytic model 
for solenoidal forcing. The analysis in Sect 14.21 showed that 
the cores tend to be smaller for compressively driven tur- 
bulence. In accordance with the clump-finding results, the 
Hennebelle-Chabrier theory implies a significant flattening 
of the high-mass wing in this case, and the slope is found to 
be close to the Salpeter value. For solenoidal forcing, how- 
ever, a significant number of cores have sizes that are of the 
same order of magnitude as the integral scale L (see Table[3] 
and Fig. |4]). This entails a violation of the basic assump- 
tion that was made in the derivation of equation (1121) . Apart 
from that, it is important to realize that both the clump- 
finding algorithm and the semi-analytic approach rely on 
the notions of the Bonnor-Ebert or Jeans mass and turbu- 
lence pressure. Being based on the collapse of a spherical 
cloud, the former is highly idealised and does not apply to 
the fully non-linear regime, while turbulence pressure is an 
ensemble-average property of isotropic turbulence. A mech- 
anism that might account for the massive cores is merging. 
Indeed, since the gas is less compressed, the cores are more 
extended and should merge more easily than in the com- 
pressible forcing case. Also they should live longer, as it 
takes more time for them to be assembled, and it is likely 
for these cores to be destroyed by violent passing shocks. 

6. Conclusions 

Computing the distributions of gravitationally unstable 
cores by means of a clump-finding algorithm as well as semi- 
analytic methods for hydrodynamic simulations of forced 
supersonic turbulence, we have found significant differences 
resulting from the turbulence forcing, the choice of the 
global mass scale, and the effects of turbulent support. For 
a comparison with theoretical predictions, the numerical 
probability density functions of the mass density fluctua- 
tions were used to evaluate the formulae for the mass dis- 
tribution. Our analysis completes a comprehensive study of 
the influence of forcing on the density an d velocity statis- 
tics of isothermal supersoii i c turbulence (iFederrath et al.l 
I2008bt iSchmidt et all l2008t IFederrath et al.l l2009aHbD ~ 
the following, we summarize the main results: 

1. For hydrodynamic simulations without explicit treat- 
ment of self-gravity, the CMF depends on the choice 
of the global mass scale, which determines the num- 
ber of Bonnor-Ebert masses, A^be, with respect to the 
mean density in the computational domain. For clump- 
finding, one has to observe two constraints: A'be must 
be sufficiently large to allow for many cores, but if 



A'be is chosen too large, a significant fraction of cores 
will be numerically unresolved and the resulting CMF 
will be shifted towards higher masses. In our analysis, 
these constraints are well satisfied for two cases: Firstly, 
solenoidal forcing with A^be ~ 1-2 x 10^ and, secondly, 
compressive forcing with A'be = 10^. The lower value 
of A^BE is consistent with Larson-type relations. 

2. Comparing the results for solenoidal and compressive 
forcing, the most noticeable trends are that the CMF 
for compressively driven turbulence displays more low- 
mass cores, while the high-mass tails are flatter, partic- 
ularly, if turbulent support is signiflcant. These trends 
are implied by equations (O and ([TU]) for the larger den- 
sity contrast produced by compressive forcing (also see 
Fig. [7]). We also found that the CMF peaks at lower 
mass in the case of compressive forcing. These results 
follow both from the clump-finding analysis and the 
semi-analytic computation of the mass distributions. 

3. Because of the considerable scatter of core masses and 
possible biases of the clump-finding algorithm, it is dif- 
ficult to ascertain power laws for the mass distributions. 
Nevertheless, we attempted to determine power-law ex- 
ponents for the high-mass tails from least square fits 
and by means of MLE. The results agree within the 
statistical uncertainties for both methods. However, the 
results fo r purely thermal support ar e at odds with the 
theory of lPadoan fc Nordlundl ()2002[ ). which asymptot- 
ically applies to A^be — >■ oo and predicts power-law tails 
for the CMF. In particular, the difference between the 
exponents following from the clump-finding analysis is 
more pronounced than the theoretical prediction. 

4. The mass distributions for cores with purely 
thermal support agree we ll with the theory of 
iHennebeUe fc Cha bricr (2008) for A^be = 10^ afthough 
the tails are less stiff towards high masses than what 
is expected on the basis of the semi-analytic models. 
There are large discrepancies for A'^be = 1-2 x 10^, 
but the disagreement might be spurious because of 
the strong dependence of the high-mass wings on the 
numerical resolution. Apparently, the mass distribu- 
tions for high values of A ^b e are in closer agreement 
with iPadoan fc Nordlundl (|2002[ ) if purely thermal 
support is assumed. However, it should be noted that 
the mass-density pdfs will be stronlgy affected by 
self-gravity in this case. 

5. For A^BE = 10'^, turbulent pressure is important for a 
wider range of core masses than in the case A^be = 1.2 x 
10^. For this reason, the CMFs change substantially 
if turbulent support is included. The effect is partic- 
ularly strong for turbulence that is driven by solenoidal 
forcing, because the mass and the size of virtually all 
cores is increased by turbulent support. However, the 
mass distribution obtained by clump-finding cannot be 
reproduced theoretically, because the basic assumption 
that the size of the cores is much smaller than the inte- 
gral scale is clearly not fulfilled in this case. Moreover, 
a stability criterion that is based on the notions of 
Jeans mass and turbulent pressure might fail if the 
cores become too large. On the other hand, we find very 
good agreement between the clump-finding analysis and 
the Hennebelle-Chabrier for compressively driven tur- 
bulence with A'be = 10^. Compared to purely thermal 
support, the high- mass tail is considerably fiatter in this 
case. 
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One of the difficulties in assessing the applicability of 
the semi-analytic models is that we cannot disentangle lim- 
itations of the underlying theoretical concepts from short- 
comings of the clump-finding algorithm. Studying the in- 
fluence of the clump-finding parameters, it becomes clear 
that the shape of the resulting mass distribution varies sig- 
nificantly with these parameters. While the peak position 
appears to be quite robust, the tails are more affected by 
the tuning of the clump-finding algorithm. Apart from these 
uncertainties and the large scatter of the core masses, there 
might be systematic biases. In particular, it is not entirely 
clear whether clump-finding traces down the smallest grav- 
itationally unstable cores. 

Even so, we have been able to shed more light on the 
gravitational fragmentation of turbulent gas. The influence 
of the forcing is irrefutable for a range of length scales that 
certainly extends beyond the energy-containing range. An 
open question is which regime in the interstellar medium is 
suitably modelled by a particular mode of forcing. If com- 
pressive excitation of turbulence occurs on length scales 
that tend to be larger than the size of molecular clouds, then 
the core statistics on much smaller scales, i. e., inside the 
cores o f the clouds, still might be universal. iFederrath et al.l 
120093), however, suggested that different forcing mecha- 
nisms can produce genuine differences in molecular clouds, 
which affect the properties of turbulence even on the length 
scales of molecular cloud cores. This calls for further ad- 
vances both on the observational and the theoretical side. 

A further important result is the influence of turbulent 
pressure on the gas fragmentation. From the theoretical 
point of view this is not new at all. But we have been able 
to demonstrate numerically that turbulent support entails 
a flattening of the CMP towards high masses. In one ex- 
treme case, the slope of the high-mass tail was found to be 
close to the Salpeter slope. In this respect, the effect of tur- 
bulent pressure on gravitational fragmentation is analogous 
to the effect of magnetic pressure in magnetohydrodynam- 
ical turbulence. This can be understood on the basis of the 
dispersion relation resulting from a linear stability analy- 
sis, which shows that the turbulent pressure as well as the 
magnetic pressure contribute to the effective pressure of 
the gas. We do not suggest that the observed CMF might 
be explainable in terms of hydrodynamical turbulence only, 
because it is known that magnetic fields play an important 
role in the interstellar medium. However, our results show 
that turbulent pressure might significantly contribute to the 
slope of the high-mass tail of the CMF. Apart from that, 
this effect is important in numerical simulations, where un- 
resolved turbulent velocity fluctuations are of the order of 
the speed of sound or higher. In this case, the corresponding 
turbule nt pressure ca n be computed with a subgrid-scale 
model (lSchmidtll2009l) . 

In order to achieve further progress with the theoretical 
explanation of the CMF, it is paramount to account for pro- 
cesses that modify the PDFs of the mass density. One such 
process is the thermal instability induced by radiative cool- 
ing in the interstellar medium. It w as already noticed by 
iPassot fc Vazauez-Semadenil (|1998[ ) that polytropic equa- 
tions of state, which are simple models for cooling, lead to 
power-law tai ls in the density PDFs. This was confirmed by 
iLi et al.l ([2002,) for three-dimensional tu r bulen ce-in-a-box 
simulations and by lAudit fc Hennebell^ (j2009D for three- 
dimensional simulations of colliding flows with a polytropic 
exponent 7 — 0.7, although a power law was not obtained 



for thermally bistable turbulence. ISeifried et al.l (j2010t) . on 
the other hand, find power-law tails for thermally bistable 
turbulence pro duced by compressive forcin g, using the cool- 
ing function of lAudit fc Hennebelli (|2005[ ). The effects of a 
polytropic equation of state o n the CMF were analysed by 
iHennebelle fc Chabried (l2009l) . 

Another process that gives rise to density PDFs 
with power-law tails is the accretion an d contraction 
of gas in sclf-gravitating turbulence (e. g. iKlessen et al.l 
2000; Federrath et al. 2008a). Moreover, dense structures 
can merge due to their mutual gravitational attraction. 
Consequentl y, the statist ical characteristics will evolve 
with time (KlessenI l2000t ). and the corresponding core 
mass distribution can differ significantly from the non-self- 
gravitating case (Klesscn 2001,), depending on the rela- 
tive strength of self-gravity and the evolutionary stage of 
the star-forming cloud. Indeed, power-law high-density tails 
are observed in h igh-resolution extinction maps of nearby 
molecular clouds (jKainulainen et aLll2009t) . 

In numerical simulations of self-gravitating turbulence, 
the mass distribution of gravitationally collapsing objects 
can be determined directly. However, the sink particles that 
are us ually app lied to numerically capture the collapsing 
gas (B ate etalirig gS: Krumholz et al. 2004 : 1 Jappsen et al.l 
120051 : iPadoan fc Nordlund.2009i : ,Federrath et al.ll2010D do 
not directly correspond to the cores we have considered 
here. While sink particles are dynamic objects that ac- 
crete gas and thus can acquire different masses, they are 
obviously not associated with a variable length scale, as 
is the case for corefl Theoretically, the relation between 
length scale and mass becomes apparent in the Hennebelle- 
Chabrier theory, but also the clump-finding algorithm iden- 
tifies regions of variable size as cores. For this reason, re- 
lating the mass-distributions of sink particles and cores is 
non-trivial. Nevertheless, the approach via sink particles 
has the considerable advantage that concepts such as the 
Bonner-Ebert mass, which is rather arbitrary in the highly 
non-linear regime, are not required. In any case, includ- 
ing self-gravity in turbulence simulations is indispensable 
to improve our understanding of gravoturbulent fragmen- 
tation. 
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^ There is a fixed accretion radius extending over a few grid 
cells, which can be interpreted as the length scale associated 
with sink particles. 
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Appendix A: Maximum likelihood estimation for 
power law scaling parameters 

Estimating the parameters a and Xmin of a power law dis- 
tribution of the form 



d{x) — Cx a; > Xn 



(A.l) 



with a constant C and exponent a by simply fitting 
straight lines on a log-log plot intr oduces large biases (e.g. 
,Newmanii2005t IClauset et al.ll2007[ ). Furthermore, this tech- 
nique gives no information whether assuming a power law 
is reasonable or not. In order to accurately estimate the 
scaling parameters of the high-mass range in our mass dis- 
tribut ions we followed the MLE method of IClauset et al.l 
(|2007[ ). In this appendix we give a brief overview of the 
method. For a detailed description please refer to the orig- 
inal paper. 

Clauset et al.l (|2007l) use a maximum likelihood estima- 
tor 



a = \ 



Vln^ 



(A.2) 



for the power law exponent a assuming that the dataset, 
consisting of n values Xi with i — 1 . . . n, is drawn from 
a continuous power law distribution for Xi > Xmin- o. is 
the estimated value for a of the underlying distribution 
(eqn. lA.ll) . The value of the lower bound Xmin of the distri- 
bution is the crucial factor in this calculation. If the esti- 
mated Xmin is too Small (xmin < Xmin) ouc will fit a powcr- 
law to non-power-law data and the parameter estimation 
will be biased. A too high imin > a:min value will throw 
away a certain amount of legitimate, power-law-distributed 
data points. A reasonable value for Xmin minimises the 
Kolmogorov-Smirnov (KS) statistics defined by 



max \S{x) 



(A.3) 



where >S'(a;) is the cumulative distribution function (CDF) 
of the observed data points and P{x) is the CDF of the 
fitted power law model. The optimal value for Xmin gives the 
minimum difference between the data and the fitted model 
which is just the minimised KS-statistics D of eqn. (|A.3p . 

To check if the observed dataset is likely to be drawn 
from a power law distribution we calculate a p- value using a 
semi-parametric bootstrap approach. Therefore, we create 
2500 synthetic datasets per snapshot which follow a real 
power law with the estimated parameters Xmin and a of 
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the observed core mass values. For xi < Xmin the synthetic 
distributions follow the same non-power law behaviour as 
the observed datasct. For each of the synthetic distributions 
we again estimate the parameters a and Xmin and compute 
the KS-statistics. The p-vahic; is the fraction of the KS- 
statistics of the synthetic datasets whose value is higher 
than the KS-statistics of the observed dataset.That means 
that p is the probability of getting a goodness of fit, e.g. 
the KS-statistics, for a real power law distributed synthetic 
dataset which is at least as bad as the goodness of fit of our 
observed dataset. For p < 0.1 the power law assumption has 
to be ruled out. Note that a value higher than 0.1 does not 
necessarily mean that the underlying distribution follows a 
power law. The p-value can only rule out the hypothesis of 
a power law distribution. 



